<?php

require_once "../Matrix.php";

class TestMatrix
{

    function TestMatrix()
    {

        // define test variables

        $errorCount = 0;
        $warningCount = 0;
        $columnwise = array(1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.);
        $rowwise = array(1., 4., 7., 10., 2., 5., 8., 11., 3., 6., 9., 12.);
        $avals = array(array(1., 4., 7., 10.), array(2., 5., 8., 11.), array(3., 6., 9., 12.));
        $rankdef = $avals;
        $tvals = array(array(1., 2., 3.), array(4., 5., 6.), array(7., 8., 9.), array(10., 11., 12.));
        $subavals = array(array(5., 8., 11.), array(6., 9., 12.));
        $rvals = array(array(1., 4., 7.), array(2., 5., 8., 11.), array(3., 6., 9., 12.));
        $pvals = array(array(1., 1., 1.), array(1., 2., 3.), array(1., 3., 6.));
        $ivals = array(array(1., 0., 0., 0.), array(0., 1., 0., 0.), array(0., 0., 1., 0.));
        $evals = array(array(0., 1., 0., 0.), array(1., 0., 2.e-7, 0.), array(0., -2.e-7, 0., 1.), array(0., 0., 1., 0.));
        $square = array(array(166., 188., 210.), array(188., 214., 240.), array(210., 240., 270.));
        $sqSolution = array(array(13.), array(15.));
        $condmat = array(array(1., 3.), array(7., 9.));
        $rows = 3;
        $cols = 4;
        $invalidID = 5; /* should trigger bad shape for construction with val        */
        $raggedr = 0; /* (raggedr,raggedc) should be out of bounds in ragged array */
        $raggedc = 4;
        $validID = 3; /* leading dimension of intended test Matrices               */
        $nonconformld = 4; /* leading dimension which is valid, but nonconforming       */
        $ib = 1; /* index ranges for sub Matrix                               */
        $ie = 2;
        $jb = 1;
        $je = 3;
        $rowindexset = array(1, 2);
        $badrowindexset = array(1, 3);
        $columnindexset = array(1, 2, 3);
        $badcolumnindexset = array(1, 2, 4);
        $columnsummax = 33.;
        $rowsummax = 30.;
        $sumofdiagonals = 15;
        $sumofsquares = 650;

        /**
         * Test matrix methods
         */

        /**
         * Constructors and constructor-like methods:
         *
         *   Matrix(double[], int)
         *   Matrix(double[][])
         *   Matrix(int, int)
         *   Matrix(int, int, double)
         *   Matrix(int, int, double[][])
         *   constructWithCopy(double[][])
         *   random(int,int)
         *   identity(int)
         */
        echo "<p>Testing constructors and constructor-like methods...</p>";

        $A = new Matrix($columnwise, 3);
        if ($A instanceof Matrix) {
            $this->try_success("Column-packed constructor...");
        } else
            $errorCount = $this->try_failure($errorCount, "Column-packed constructor...", "Unable to construct Matrix");

        $T = new Matrix($tvals);
        if ($T instanceof Matrix)
            $this->try_success("2D array constructor...");
        else
            $errorCount = $this->try_failure($errorCount, "2D array constructor...", "Unable to construct Matrix");

        $A = new Matrix($columnwise, $validID);
        $B = new Matrix($avals);
        $tmp = $B->get(0, 0);
        $avals[0][0] = 0.0;
        $C = $B->minus($A);
        $avals[0][0] = $tmp;
        $B = Matrix::constructWithCopy($avals);
        $tmp = $B->get(0, 0);
        $avals[0][0] = 0.0;
        /** check that constructWithCopy behaves properly **/
        if (($tmp - $B->get(0, 0)) != 0.0)
            $errorCount = $this->try_failure($errorCount, "constructWithCopy... ", "copy not effected... data visible outside");
        else
            $this->try_success("constructWithCopy... ", "");

        $I = new Matrix($ivals);
        if ($this->checkMatrices($I, Matrix::identity(3, 4)))
            $this->try_success("identity... ", "");
        else
            $errorCount = $this->try_failure($errorCount, "identity... ", "identity Matrix not successfully created");

        /**
         * Access Methods:
         *
         *   getColumnDimension()
         *   getRowDimension()
         *   getArray()
         *   getArrayCopy()
         *   getColumnPackedCopy()
         *   getRowPackedCopy()
         *   get(int,int)
         *   getMatrix(int,int,int,int)
         *   getMatrix(int,int,int[])
         *   getMatrix(int[],int,int)
         *   getMatrix(int[],int[])
         *   set(int,int,double)
         *   setMatrix(int,int,int,int,Matrix)
         *   setMatrix(int,int,int[],Matrix)
         *   setMatrix(int[],int,int,Matrix)
         *   setMatrix(int[],int[],Matrix)
         */
        print "<p>Testing access methods...</p>";

        $B = new Matrix($avals);
        if ($B->getRowDimension() == $rows)
            $this->try_success("getRowDimension...");
        else
            $errorCount = $this->try_failure($errorCount, "getRowDimension...");

        if ($B->getColumnDimension() == $cols)
            $this->try_success("getColumnDimension...");
        else
            $errorCount = $this->try_failure($errorCount, "getColumnDimension...");

        $barray = $B->getArray();
        if ($this->checkArrays($barray, $avals))
            $this->try_success("getArray...");
        else
            $errorCount = $this->try_failure($errorCount, "getArray...");

        $bpacked = $B->getColumnPackedCopy();
        if ($this->checkArrays($bpacked, $columnwise))
            $this->try_success("getColumnPackedCopy...");
        else
            $errorCount = $this->try_failure($errorCount, "getColumnPackedCopy...");

        $bpacked = $B->getRowPackedCopy();
        if ($this->checkArrays($bpacked, $rowwise))
            $this->try_success("getRowPackedCopy...");
        else
            $errorCount = $this->try_failure($errorCount, "getRowPackedCopy...");

        /**
         * Array-like methods:
         *   minus
         *   minusEquals
         *   plus
         *   plusEquals
         *   arrayLeftDivide
         *   arrayLeftDivideEquals
         *   arrayRightDivide
         *   arrayRightDivideEquals
         *   arrayTimes
         *   arrayTimesEquals
         *   uminus
         */
        print "<p>Testing array-like methods...</p>";

        /**
         * I/O methods:
         *   read
         *   print
         *   serializable:
         *   writeObject
         *   readObject
         */
        print "<p>Testing I/O methods...</p>";

        /**
         * Test linear algebra methods
         */
        echo "<p>Testing linear algebra methods...<p>";

        $A = new Matrix($columnwise, 3);
        if ($this->checkMatrices($A->transpose(), $T))
            $this->try_success("Transpose check...");
        else
            $errorCount = $this->try_failure($errorCount, "Transpose check...", "Matrices are not equal");

        if ($this->checkScalars($A->norm1(), $columnsummax))
            $this->try_success("Maximum column sum...");
        else
            $errorCount = $this->try_failure($errorCount, "Maximum column sum...", "Incorrect: " . $A->norm1() . " != " . $columnsummax);

        if ($this->checkScalars($A->normInf(), $rowsummax))
            $this->try_success("Maximum row sum...");
        else
            $errorCount = $this->try_failure($errorCount, "Maximum row sum...", "Incorrect: " . $A->normInf() . " != " . $rowsummax);

        if ($this->checkScalars($A->normF(), sqrt($sumofsquares)))
            $this->try_success("Frobenius norm...");
        else
            $errorCount = $this->try_failure($errorCount, "Frobenius norm...", "Incorrect:" . $A->normF() . " != " . sqrt($sumofsquares));

        if ($this->checkScalars($A->trace(), $sumofdiagonals))
            $this->try_success("Matrix trace...");
        else
            $errorCount = $this->try_failure($errorCount, "Matrix trace...", "Incorrect: " . $A->trace() . " != " . $sumofdiagonals);

        $B = $A->getMatrix(0, $A->getRowDimension(), 0, $A->getRowDimension());
        if ($B->det() == 0)
            $this->try_success("Matrix determinant...");
        else
            $errorCount = $this->try_failure($errorCount, "Matrix determinant...", "Incorrect: " . $B->det() . " != " . 0);

        $A = new Matrix($columnwise, 3);
        $SQ = new Matrix($square);
        if ($this->checkMatrices($SQ, $A->times($A->transpose())))
            $this->try_success("times(Matrix)...");
        else {
            $errorCount = $this->try_failure($errorCount, "times(Matrix)...", "Unable to multiply matrices");
            $SQ->toHTML();
            $AT->toHTML();
        }

        $A = new Matrix($columnwise, 4);

        $QR = $A->qr();
        $R = $QR->getR();
        $Q = $QR->getQ();
        if ($this->checkMatrices($A, $Q->times($R)))
            $this->try_success("QRDecomposition...", "");
        else
            $errorCount = $this->try_failure($errorCount, "QRDecomposition...", "incorrect qr decomposition calculation");

        $A = new Matrix($columnwise, 4);
        $SVD = $A->svd();
        $U = $SVD->getU();
        $S = $SVD->getS();
        $V = $SVD->getV();
        if ($this->checkMatrices($A, $U->times($S->times($V->transpose()))))
            $this->try_success("SingularValueDecomposition...", "");
        else
            $errorCount = $this->try_failure($errorCount, "SingularValueDecomposition...", "incorrect singular value decomposition calculation");

        $n = $A->getColumnDimension();
        $A = $A->getMatrix(0, $n - 1, 0, $n - 1);
        $A->set(0, 0, 0.);

        $LU = $A->lu();
        $L = $LU->getL();
        if ($this->checkMatrices($A->getMatrix($LU->getPivot(), 0, $n - 1), $L->times($LU->getU())))
            $this->try_success("LUDecomposition...", "");
        else
            $errorCount = $this->try_failure($errorCount, "LUDecomposition...", "incorrect LU decomposition calculation");

        $X = $A->inverse();
        if ($this->checkMatrices($A->times($X), Matrix::identity(3, 3)))
            $this->try_success("inverse()...", "");
        else
            $errorCount = $this->try_failure($errorCount, "inverse()...", "incorrect inverse calculation");

        $DEF = new Matrix($rankdef);
        if ($this->checkScalars($DEF->rank(), min($DEF->getRowDimension(), $DEF->getColumnDimension()) - 1))
            $this->try_success("Rank...");
        else
            $this->try_failure("Rank...", "incorrect rank calculation");

        $B = new Matrix($condmat);
        $SVD = $B->svd();
        $singularvalues = $SVD->getSingularValues();
        if ($this->checkScalars($B->cond(), $singularvalues[0] / $singularvalues[min($B->getRowDimension(), $B->getColumnDimension()) - 1]))
            $this->try_success("Condition number...");
        else
            $this->try_failure("Condition number...", "incorrect condition number calculation");

        $SUB = new Matrix($subavals);
        $O = new Matrix($SUB->getRowDimension(), 1, 1.0);
        $SOL = new Matrix($sqSolution);
        $SQ = $SUB->getMatrix(0, $SUB->getRowDimension() - 1, 0, $SUB->getRowDimension() - 1);
        if ($this->checkMatrices($SQ->solve($SOL), $O))
            $this->try_success("solve()...", "");
        else
            $errorCount = $this->try_failure($errorCount, "solve()...", "incorrect lu solve calculation");

        $A = new Matrix($pvals);
        $Chol = $A->chol();
        $L = $Chol->getL();
        if ($this->checkMatrices($A, $L->times($L->transpose())))
            $this->try_success("CholeskyDecomposition...", "");
        else
            $errorCount = $this->try_failure($errorCount, "CholeskyDecomposition...", "incorrect Cholesky decomposition calculation");

        $X = $Chol->solve(Matrix::identity(3, 3));
        if ($this->checkMatrices($A->times($X), Matrix::identity(3, 3)))
            $this->try_success("CholeskyDecomposition solve()...", "");
        else
            $errorCount = $this->try_failure($errorCount, "CholeskyDecomposition solve()...", "incorrect Choleskydecomposition solve calculation");

        $Eig = $A->eig();
        $D = $Eig->getD();
        $V = $Eig->getV();
        if ($this->checkMatrices($A->times($V), $V->times($D)))
            $this->try_success("EigenvalueDecomposition (symmetric)...", "");
        else
            $errorCount = $this->try_failure($errorCount, "EigenvalueDecomposition (symmetric)...", "incorrect symmetric Eigenvalue decomposition calculation");

        $A = new Matrix($evals);
        $Eig = $A->eig();
        $D = $Eig->getD();
        $V = $Eig->getV();
        if ($this->checkMatrices($A->times($V), $V->times($D)))
            $this->try_success("EigenvalueDecomposition (nonsymmetric)...", "");
        else
            $errorCount = $this->try_failure($errorCount, "EigenvalueDecomposition (nonsymmetric)...", "incorrect nonsymmetric Eigenvalue decomposition calculation");

        print("<b>{$errorCount} total errors</b>.");
    }

    /**
     * Print appropriate messages for successful outcome try
     * @param string $s
     * @param string $e
     */
    function try_success($s, $e = "")
    {
        print "> " . $s . "success<br />";
        if ($e != "")
            print "> Message: " . $e . "<br />";
    }

    /**
     * Print appropriate messages for unsuccessful outcome try
     * @param int $count
     * @param string $s
     * @param string $e
     * @return int incremented counter
     */
    function try_failure($count, $s, $e = "")
    {
        print "> " . $s . "*** failure ***<br />> Message: " . $e . "<br />";
        return ++$count;
    }

    /**
     * Print appropriate messages for unsuccessful outcome try
     * @param int $count
     * @param string $s
     * @param string $e
     * @return int incremented counter
     */
    function try_warning($count, $s, $e = "")
    {
        print "> " . $s . "*** warning ***<br />> Message: " . $e . "<br />";
        return ++$count;
    }

    /**
     * Check magnitude of difference of "scalars".
     * @param float $x
     * @param float $y
     */
    function checkScalars($x, $y)
    {
        $eps = pow(2.0, -52.0);
        if ($x == 0 & abs($y) < 10 * $eps) return;
        if ($y == 0 & abs($x) < 10 * $eps) return;
        if (abs($x - $y) > 10 * $eps * max(abs($x), abs($y)))
            return false;
        else
            return true;
    }

    /**
     * Check norm of difference of "vectors".
     * @param float $x []
     * @param float $y []
     */
    function checkVectors($x, $y)
    {
        $nx = count($x);
        $ny = count($y);
        if ($nx == $ny)
            for ($i = 0; $i < $nx; ++$i)
                $this->checkScalars($x[$i], $y[$i]);
        else
            die("Attempt to compare vectors of different lengths");
    }

    /**
     * Check norm of difference of "arrays".
     * @param float $x [][]
     * @param float $y [][]
     */
    function checkArrays($x, $y)
    {
        $A = new Matrix($x);
        $B = new Matrix($y);
        return $this->checkMatrices($A, $B);
    }

    /**
     * Check norm of difference of "matrices".
     * @param matrix $X
     * @param matrix $Y
     */
    function checkMatrices($X = null, $Y = null)
    {
        if ($X == null || $Y == null)
            return false;

        $eps = pow(2.0, -52.0);
        if ($X->norm1() == 0. & $Y->norm1() < 10 * $eps) return true;
        if ($Y->norm1() == 0. & $X->norm1() < 10 * $eps) return true;

        $A = $X->minus($Y);

        if ($A->norm1() > 1000 * $eps * max($X->norm1(), $Y->norm1()))
            die("The norm of (X-Y) is too large: " . $A->norm1());
        else
            return true;
    }

}

$test = new TestMatrix;
?>
